function [angle0_output, angle6_result, compound_result, comparison_result] = CompoundPlaneWaveSimulation()
% %输入参数：
% %无
% %输出参数：
% % angle0_output      ：只使用第0角度的成像结果，图像结构体，即图像及图像相关参数mat
% % angle6_result      ：只使用第6角度的成像结果，图像结构体，即图像及图像相关参数mat
% % compound_result    : 多个角度平面波复合成像的结果，图像结构体，即图像及图像相关参数mat
% % comparison_result  : 单角度平面波成像结果与多角度平面波成像结果的比较，两条曲线

%=============Define the tubes in phantom=============%
xs = [-2.1:0.1:-1.5 -0.9:0.1:-0.3 0.3:0.1:0.9 1.5:0.1:2.1]*1e-3; % in mm
xs = repmat(xs,6,1); % in mm
zs = [5.6;5.7;5.8;5.9;6.0;6.1]*1e-3;% in mm
zs = repmat(zs,1,28);
RC = ones(size(xs));
%=============Define the tubes in phantom=============%

%=============Define Transducer param=================%
param.fc = 2.8*1e6; %[Hz]
param.kerf = 0.046*1e-3; % kerf == [mm]   
param.width = 0.254*1e-3;  % widch == [mm]   
param.pitch = 0.3*1e-3; % pitch == [mm]  pitch need to equal to width + kerf 
param.Nelements = 64;
param.bandwidth = 90; % [percent]
param.radius = Inf;
param.fs = 4*param.fc; % sampling frequency i Hz
param.c = 1500; % [m/s], the reconstruction sound speed - water 1500 m/s
param.fnumber = []; % determine automatically

%=============Define Transducer param=================%

%=============Define plane-wave======================%

titlangle = 15;  % tiltangles value {\circ}
na=11; % Define the numbers of tilt angles, it must be an odd
tilt = linspace(-titlangle,titlangle,na)/180*pi; % -tiltangle 0 +tiltangle
txdelayindex = cell(na,1); % this cell will contain the transmit delays index
resoltionpamater = 300;    % typecal 200
log_compressed_parameter = 36; % log-compressed-parameter dB typecal 36
%=============Define plane-wave======================%

%=============Generation plane-wave and sampling======================%
for k = 1:na
    txdelayindex{k} = txdelay(param,tilt(k));
end

RF = cell(na,1); % this cell will contain the RF with na angles
option.WaitBar = false; % remove the progress bar of SIMUS
for k = 1:na
    RF{k} = simus(xs,zs,RC,txdelayindex{k},param,option);
end

IQ = cell(na,1);  % this cell will contain the I/Q

for k = 1:na
    IQ{k} = rf2iq(RF{k},param.fs,param.fc);
end

param.fnumber = [];
[xi,zi] = meshgrid(linspace(-6e-3,6e-3,resoltionpamater),linspace(0,12e-3,resoltionpamater));% spaceing = 12e-3/300
beamformed_IQ = zeros(resoltionpamater,resoltionpamater,na); % this array will contain the na I/Q images

%h = waitbar(0,'');
for k = 1:na
%    waitbar(k/na,h,['DAS processing: I/Q # ' int2str(k) ' of total angles # '  int2str(na) ])
    beamformed_IQ(:,:,k) = das(IQ{k},xi,zi,txdelayindex{k},param);
end
%close(h)

%% angle 0

angle0_output.image = bmode(beamformed_IQ(:,:,1),log_compressed_parameter); % log-compressed image
angle0_output.image = transpose(angle0_output.image);
angle0_output.dimension = size(angle0_output.image);
angle0_output.pixel_spacing = [0.04,0.04];
angle0_output.Xposition = xi(1,:)*1e3; %imagesc(angle0_output.Xposition,angle0_output.Zposition,angle0_output.image);
angle0_output.Zposition = zi(:,1)*1e3;
angle0_output.datatype = class(angle0_output.image);
%%angle0_output.Colormap = get(gca,'Colormap');
angle0_output.XTickLabel = '[mm]';
angle0_output.YTickLabel = '[mm]';
angle0_output.title = 'DAS image with the 0st tilt angle';

%%imwrite(angle0_output.image, 'output.bmp');

close(gcf);
%% angle 6

angle6_result.image = bmode(beamformed_IQ(:,:,(na-1)/2+1),log_compressed_parameter); % log-compressed image
angle6_result.image = transpose(angle6_result.image);
angle6_result.dimension = size(angle6_result.image);
angle6_result.pixel_spacing = [0.04,0.04];
angle6_result.Xposition = xi(1,:)*1e3;
angle6_result.Zposition = zi(:,1)*1e3;
angle6_result.datatype = class(angle6_result.image);
%%angle6_result.Colormap = get(gca,'Colormap');
angle6_result.XTickLabel = '[mm]';
angle6_result.YTickLabel = '[mm]';
angle6_result.title = 'DAS image with the 6th tilt angle';
close(gcf);
%% compound_result

log_compressed_parameter=8;
cIQ = sum(beamformed_IQ,3); % Compound all angles beamformed I/Q

compound_result.image = bmode(cIQ,log_compressed_parameter); % log-compressed image
compound_result.image = transpose(compound_result.image);
compound_result.dimension = size(compound_result.image);
compound_result.pixel_spacing = [0.04,0.04];
compound_result.Xposition = xi(1,:)*1e3;
compound_result.Zposition = zi(:,1)*1e3;
compound_result.datatype = class(compound_result.image);
%%compound_result.Colormap = get(gca,'Colormap');
compound_result.XTickLabel = '[mm]';
compound_result.YTickLabel = '[mm]';
compound_result.title = 'PW-based DAS image with compound 11 angles';
close(gcf);

%% comparison_result

line1 = bmode(beamformed_IQ(:,:,(na-1)/2+1),log_compressed_parameter);
line2 = bmode(cIQ,log_compressed_parameter);

comparison_result.OneAngle = line1(150,:);
comparison_result.compound = line2(150,:);
comparison_result.XAxisTickLabel = [0 50 100 150 200 250 300];
comparison_result.YAxisTickLabel = [0 50 100 150];
comparison_result.XAxisLabel = 'time (μs)';
comparison_result.YAxisLabel = 'Amplitude (a.u.)';
comparison_result.legend = {'DAS with only one angle','PW-based DAS image'};
comparison_result.title = 'RF from 64th element (0st tilt angle, 1st frame)';

end





% function [angle1_result, angle6_result, compound_result, comparison_result] = CompoundPlaneWaveSimulation()
% %输入参数：
% %无
% %输出参数：
% % angle0_output      ：只使用第0角度的成像结果，图像结构体，即图像及图像相关参数mat
% % angle6_result      ：只使用第6角度的成像结果，图像结构体，即图像及图像相关参数mat
% % compound_result    : 多个角度平面波复合成像的结果，图像结构体，即图像及图像相关参数mat
% % comparison_result  : 单角度平面波成像结果与多角度平面波成像结果的比较，两条曲线
% 
% %% angle1_output, 一个图像（mat）需包括如下元素（image, dimension, pixel_spacing, channels, datatype, Colormap, XTick, XTickLabel, YTick, YTickLabel）用于前端显示
% angle1_img = imread("../const_file/one_angle_first.bmp");
% angle1_result.image = angle1_img;
% angle1_result.dimension = size(angle1_img);
% angle1_result.pixel_spacing = [0.05,0.05];
% angle1_result.channels = 3;
% angle1_result.datatype = class(angle1_img);
% angle1_result.Colormap = transpose(colormap("hot"));
% angle1_result.XTick = [1 96 192 288 384];
% angle1_result.XTickLabel = '[mm]';
% angle1_result.YTick = [1 100 200 300 400 500 600 700 800 900];
% angle1_result.YTickLabel = '[mm]';
% %% 如下为mock结果
% fid = fopen("../const_file/one_angle_first.bmp", "r", "b");
% file_content = fread(fid);
% angle1_result.image = cast(file_content, "uint8");
% fclose(fid);
% 
% %%  angle6_result, 一个图像（mat）需包括如下元素用于前端显示
% angle6_img = imread("../const_file/one_angle_sixth.bmp");
% angle6_result.image = angle6_img;
% angle6_result.dimension = size(angle6_img);
% angle6_result.pixel_spacing = [0.05,0.05];
% angle6_result.channels = 3;
% angle6_result.datatype = class(angle6_img);
% angle6_result.Colormap = transpose(colormap("hot"));
% angle6_result.XTick = [1 96 192 288 384];
% angle6_result.XTickLabel = '[mm]';
% angle6_result.YTick = [1 100 200 300 400 500 600 700 800 900];
% angle6_result.YTickLabel = '[mm]';
% %% 如下为mock结果
% fid = fopen("../const_file/one_angle_sixth.bmp", "r", "b");
% file_content = fread(fid);
% angle6_result.image = cast(file_content, "uint8");
% fclose(fid);
% 
% %%  compound_result, 一个图像（mat）需包括如下元素用于前端显示
% compound_img = imread("../const_file/compound_plane_wave.bmp");
% %imshow(img)
% compound_result.image = compound_img;
% compound_result.dimension = size(compound_img);
% compound_result.pixel_spacing = [0.05,0.05];
% compound_result.channels = 3;
% compound_result.datatype = class(compound_img);
% compound_result.Colormap = transpose(colormap("hot"));
% compound_result.XTick = [1 96 192 288 384];
% compound_result.XTickLabel = '[mm]';
% compound_result.YTick = [1 100 200 300 400 500 600 700 800 900];
% compound_result.YTickLabel = '[mm]';
% %% 如下为mock结果
% fid = fopen("../const_file/compound_plane_wave.bmp", "r", "b");
% file_content = fread(fid);
% compound_result.image = cast(file_content, "uint8");
% fclose(fid);
% 
% %%  comparison_result, 单角度平面波成像结果与多角度平面波成像结果的比较，两条曲线
% x = 0:0.01:3;
% comparison_result.one_angle_signals = sin(x);
% comparison_result.one_angle_signals_title = "DAS with only one angle";
% comparison_result.compound_signals = cos(x);
% comparison_result.compound_signals_title = "PW-based DAS image";
% comparison_result.XAxisTickLabel = [0 50 100 100 150 200 250 300];
% comparison_result.XAxisLabel = '[time(us)]';
% comparison_result.YAxisTickLabel = [0 50 100 150 900];
% comparison_result.YAxisLabel = '[Amplitude(a.u.)]';
% 
% end

